Transit facility allocation: Hybrid quantum-classical optimization

An essential consideration in urban transit facility planning is service efficiency and accessibility. Previous research has shown that reducing the number of facilities along a route may increase efficiency but decrease accessibility. Striking a balance between these two is a critical consideration in transit planning. Transit facility consolidation is a cost-effective way to improve the quality of service by strategically determining the desirable allocation of a limited number of facilities. This paper develops an optimization framework that integrates Geographical Information systems (GIS), decision-making analysis, and quantum technologies for addressing the problem of facility consolidation. Our proposed framework includes a novel mathematical model that captures non-linear interactions between facilities and surrounding demand nodes, inter-facility competition, ridership demand and spatial coverage. The developed model can harness the power of quantum effects such as superposition and quantum tunnelling and enables transportation planners to utilize the most recent hardware solutions such as quantum and digital annealers, coherent Ising Machines and gate-based universal quantum computers. This study presents a real-world application of the framework to the public transit facility redundancy problem in the British Columbia Vancouver metropolitan area. We demonstrate the effectiveness of our framework by reducing the number of facilities by 40% while maintaining the same service accessibility. Additionally, we showcase the ability of the proposed mathematical model to take advantage of quantum annealing and classical optimization techniques.


Introduction
In a growing society, public transportation systems are widely used to serve large populations in a space-time-concentrated travel demand manner [1]. The more successful a public transportation system is, the greater the positive effect on the population who uses it [2]. The success of a transportation system requires a balance between the accessibility and efficiency of the system [3]. Historically, some transit systems have issues achieving this balance.
Accessibility refers to how easily demand nodes (e.g. urban areas) can access the transit facilities [4]. Generally, the closer a facility is to a demand node, the more accessible the facility becomes. Thus, more facilities in place along a given route would create more accessibility for transit riders [4]. Efficiency is defined by how much distance a transit route can cover in a a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 particular time window [4]. The efficiency of a transit route depends on a few factors such as the number of facilities on the route, traffic flow, dwell times at facilities, and the number of passengers boarding and alighting at each transit facility. Ideally, a transit route is highly efficient to appease customers and costs in travel time. However, finding a balance between accessibility and efficiency is preferred. As a case study, we consider the transportation system of one of the densest cities in British Columbia, Vancouver. Translink is a mass transportation system responsible for the regional transportation network of Metro Vancouver, including public transport, major roads and bridges. We use publicly available GIS data from Translink and Census program to optimize one of the most problematic bus routes. We demonstrate that the proposed framework can reduce the number of transit facilities by 40% while maintaining the same accessibility and improving efficiency.

Contributions
With these definitions for accessibility and efficiency, the main goal is to determine the best allocation of a limited number of transit facilities. This paper proposes a novel end-to-end optimization framework that accomplishes the goal. Specifically, the paper's contributions are: • We propose a novel mathematical model inspired by the non-linear Spatial Interaction Coverage (SIC) models [5,6]. The proposed model addresses transit facility allocation and accessibility by modelling population ridership (demand), physical accessibility (distances between demand nodes and facilities), connectivity (number of destinations from a facility), distance decay and inter-facility competition. One of the model's unique features is that its solutions could be naturally integral even if the integrality constraints are disregarded. See Section 2 and Section 4.
• Although expressive, SIC models are extremely hard to solve [5,6]. We address this issue by using Quadratic Unconstrained Binary Optimisation (QUBO), which is closely related to the Ising spin glass model used in Statistical Physics and Quantum Simulation [7,8]. Due to this relation, our model is readily solvable on a wide variety of classical and quantum optimization metaheuristics such as Quantum Annealing, quantum-inspired Digital Annealing [9], quantum optics-based hardware such as Coherent Ising Machines [10] and universal gate quantum computers [11]. Thus, the proposed model enables transportation planners to utilize and experiment with quantum-based technology along with a conventional classicalcomputing toolset. See Section 3.
• To analyze the quality of attained solutions, we present a theoretical upper bound on optimal values of the proposed model. The upper bound is fairly tight and can be computed in polynomial time. See Section 5.
• We propose a flexible integration of GIS into the model using the multi-criteria decisionmaking analysis [12]. The analysis is widely used in many applications related to decisionmaking in logistics, industry and business. See Section 3.3.
In addition to the above, we present a hybrid quantum-classical solver built with the opensource dwave-hybrid Python framework [13]. To explore solution search space robustly, the solver utilizes classical hill-climbing mechanics together with quantum superposition and quantum tunnelling effects [14].
The remainder of this paper proceeds as follows. Section 1.2 reviews various optimization models and algorithms. Section 2 describes the formulation of the proposed mathematical model. Section 3 describes QUBO formalism, model simplification, application of TOPSIS and quantum-classical solver. Sections 4 and Section 5 discuss the model's properties, such as solutions' integrality and an objective function bound. Section 6 showcases the application of the framework to the transit bus route in Vancouver, British Columbia. Section 8 summarizes both mathematical and numerical findings and highlights the implications of this study. Section 9 contains concluding remarks and future research directions.

Related literature
There are several modelling approaches for locating optimal facilities. However, these models do not simultaneously address distance decay, coverage range, competitors and GIS data. A well-known approach for addressing a service coverage range is Location Set Covering Problem (LSCP) [15]. LSCP seeks to place the minimum number of facilities such that all demand nodes can be served. While LSCP can identify the minimal subset of facilities necessary for covering all demand nodes it is unable to account for distance attenuation [5]. In some application scenarios, the model's requirement for all demands to be covered can be a fairly stringent condition that limits the choice of possible facility configurations. Moreover, often we would like to control the number of chosen facilities and not necessarily aim for the smallest number of them. Another modelling method is a distance-based approach such as the pmedian problem [16]. The p-median formulation aims to minimize the total demandweighted distance between each demand node and the nearest facility. While the p-median problem implements weighted distance as one of the decision factors and allows to control the number of chosen facilities, the binary nature of facility assignment does not allow for a demand node to be served by an arbitrary number of facilities.
Furthermore, recent studies actively integrate GIS into transit optimization [17]. GIS introduces more refined and data-oriented modelling, allowing us to account for census geographies such as population density per area, dissemination area boundaries and global positioning data. One of the significant requirements of GIS approach is large-scale data handling. Studies that incorporate GIS often use closed-source and inaccessible software which imposes strict requirements on how data should be handled. This in turn makes the proposed methodology not applicable to a wider range of scenarios. One such GIS-based approach is introduced in [3,5]. The study combines the aspects of the p-median problem with the Maximal Covering Location Problem (MCLP) discussed in [5]. The proposed Spatial Interaction Coverage (SIC) model introduces the concept of interaction between a demand node and a facility. The interaction incorporates GIS-based demand, distance decay, geographical coverage and inter-facility competition. The model aims to maximize the interaction between a demand node and a facility. The resulting formulation is a non-linear, binary integer problem. The downside of the approach is the necessity to develop a proprietary optimization heuristic that needs tight coupling with commercial GIS software. The non-linearity of the problem combined with integrality constraints significantly reduces the choices of available solvers.
Several studies investigate an integrated Multiple Criteria Decision Making (MCDM) approach to optimize the facility-location problem. Such decisions require an algorithm that optimizes the quantitative and qualitative factors of each facility. These take the form of a scoring/ranking algorithm that takes into consideration factors including: a transit facility's connectedness to adjoining routes, ridership demand and its location [18]. For a comprehensive multi-criteria decision analysis [18] proposes Analytical Hierarchical Process (AHP) as such a method. AHP is a structured technique for organizing and analyzing complex decisions based on pairwise comparisons of desired criteria. AHP has several drawbacks: judgment criteria are not guaranteed to be always consistent and it does not scale well to large datasets [19]. An alternative to AHP is the Technique for Order Preference by Similarity to Ideal Solution (TOPSIS) [12,20]. TOPSIS is based on an aggregating function representing "closeness to the ideal", which originated in the compromise programming method. The idea behind TOPSIS is computing a synthetic ideal candidate and determining the closeness of each facility to the ideal. Unlike AHP, TOPSIS scales well to large datasets, it can process hundreds of candidates and dozens of decision criteria, finding an optimal compromise in difficult decision situations.
Several well-known options exist for implementing and solving combinatorial problems. Commercial solvers, such as Lingo, Gurobi or CPLEX can be used, but also programming languages like Python allow the development of custom solvers that can tackle arbitrary combinatorial problems. Using objected-oriented programming has made it easier to integrate GIS data into a model formulation [3]. As mentioned previously, our paper discusses the formulation of a mathematical model using QUBO as a framework. In that regard, there exists a wide variety of powerful software-and hardware-based QUBO solvers which are successfully applied to large-scale optimization problems. QUBO problems can be efficiently solved on classical central processing units (CPU), graphics processing units (GPU) as well as specialized optimization hardware such as D-Wave's Quantum Annealer (DA) [21,22], Fujitsu's Digital Annealer (DA) [9] and the Coherent Ising Machine (CIM) [10]. In addition to hardware solutions, there exist a number of software-based optimization algorithms that can efficiently solve QUBO problems. One of such algorithms is Tabu Search (TS). TS is a metaheuristic search method. Unlike local search methods which often get trapped in suboptimal regions, TS avoids this behaviour by prohibiting already visited solutions [23]. Simulated annealing (SA) is another probabilistic hill-climbing heuristic that searches for a global optimum [24]. Combining software and specialized hardware solutions into a single algorithm yields a powerful and versatile approach for solving optimization problems. We summarize both quantum and classical heuristics that can readily solve our proposed model in Table 1.

Model formulation
We base our formulation on the aforementioned MCLP and SIC models. To reflect the current demand along the route, we extensively utilize GIS Census Program data and work on the scale of the Dissemination Area (also known as the Census Block or Dissemination Block), which is one of the minimal geographic units in the population Census Program in North America. The objective of the facility allocation optimization problem is to select p optimal facilities such that these p facilities maximize route accessibility while reducing redundancy. We formulate the model as a maximization of a QUBO problem. All QUBO problems consist of linear and quadratic parts. As we will see later the linear part of our model accounts for relations between facilities and neighbouring demand nodes (Dissemination Areas) within a radius R 0 . And the quadratic part accounts for inter-facility competition, i.e. if a facility has neighbouring competitors within a defined radius R 1 the facility has less importance during an optimization process. This approach helps to determine optimally spaced facilities that cover the maximum possible demand along the route. Our mathematical notation is summarized in Table 2.
We commence by introducing the model formulation in a quadratic form subject to the constraint X j2F where Our goal is to maximize the quadratic function (1) by activating an optimal combination of the x j for j = 1, 2, . . ., |F| such that there are a total of p active facilities, i.e., the p-constraint (2) is satisfied. In other words, we want to select the best possible combination of p facilities that maximize the objective. To develop an intuition behind (1) we look at its first linear part ∑ j2F Q jj x j . LetQ  This quantity represents the weighted sum of population a k attenuated by distance d jk between the facility j and the demand node k. We note that a demand node k is in the neighbourhood A j of a facility j. Thus each facility j has an associated aggregated population, see Fig 1. The greater the weighted population the higher the likelihood that the facility will be selected as a part of the final solution. We now explain the meaning of (3). As mentioned previously, our model accounts for inter-facility competition. We set the following simple rule: if a facility j has one neighbouring competitor then the aggregated population can be distributed among two facilities. Thus we divideQ jj by 2. Similarly, if a facility j has 2 competitors, the aggregated population can be distributed between 3 facilities, hence we divideQ jj by 3 and so on. This relation can be represented as follows: where F j is a neighbourhood of competing facilities. For example, if all the x i are active, the aggregated population of a facility j is divided by is non-linear and hard to solve. Therefore, we approximate it via a quadratic form by Finally, we observe that this gives the coefficients of the objective function This shows that all factors involving m j come from linearization around m j of (6) and that the quadratic part of the objective function models competition among the facilities. Such a formulation has a number of interesting properties: • The optimization problem can be intuitively studied from a graph-theoretical perspective (Section 3.2).
• The formulation can be easily cast into QUBO formalism.
• Both discrete and continuous optimization methods yield solutions to QUBO that are naturally binary, i.e. under certain conditions we can disregard integrality requirements and still obtain a binary solution (Section 4).
• We can derive a fairly tight upper bound for a maximum value of the objective function (Section 5).

Small example
We now look at a small scale example with 3 facilities, i.e. F = {1, 2, 3}. Without loss of generality, consider the second facility, with j = 2. Suppose that the facilities j = 1 and j = 3 are competitors of the second facility, thus the neighbourhood And the objective function isf For instance, deactivating the second facility (x 2 = 0) completely removes contributions of a facility j = 2 to the objective function. That is,f 2 ¼ 0, while bothf 1 andf 3 have a greater contribution to the objective function value because Q 21 x 2 x 1 , Q 23 x 2 x 3 = 0. Whereas, deactivating x 1 will increase the importance of the facility j = 2. It is important to note that the competition relation is reflexive, i.e. if the facility j = 2 has the competitor j = 3 then F 3 also contains the facility j = 2. Such reflexive relation leads to a combinatorial explosion of possible configurations of x j .

Optimization methods
This section describes QUBO formalism and how to use it to obtain the Ising equivalent model. Furthermore, we demonstrate simplification of the p-constraint that reduces the complexity of the problem and provides more controlled optimization results. We conclude the section with the GIS integration using TOPSIS and the presentation of the Hybrid quantumclassical solver.

QUBO formalism
So far we have presented a quadratic objective function with one explicit p-constraint (2). In order to obtain the final QUBO formulation, we need to incorporate the constraint into the objective. Specifically, for a positive scalar γ, we subtract a quadratic penalty from the objective function. Hence, we now consider the problem We note that the p-constraint (2) is now a soft constraint, i.e. it is incorporated into the objective function as a quadratic penalty. It can be shown that (16) can be compactly expressed in a matrix form. Noting that x 2 i ¼ x i for i = 1, . . ., |F|, the quadratic objective can be rewritten in a matrix form as X where x 2 {0, 1} |F| and the matrix Q 2 R jFj�jFj has entries Q ij defined in (3) and (4). Similarly, the penalty term can be expressed as We can drop the constant term p t p = p 2 , and noting that x 2 i ¼ x i we can rewrite the above as Let G ¼ 11 t À 2pI. Then the final QUBO formulation is where � Q≔Q À gG is a symmetric matrix. Therefore, our objective is to solve the following optimization problem

A graph theoretical view and partitioning the constraint
One of the advantages of the QUBO formulation is that we can view the entire problem as a network graph and assess its complexity in a visual manner. To this end, we view each nonzero Q ij of the quadratic objective presented in (1) as an edge of a graph with |F| vertices. Since (1) relates pairs of facilities by their neighbourhood-dependent competition, we expect the resulting graph to be sparse: namely, each vertex (facility) is connected to a neighbouring vertex on a transit route. See Fig 2, left. Such a graph is almost a tree and it is easy to work with. However, upon the addition of the p-constraint (2) the resulting graph is a complete graphsee Fig 2, center. This can be intuitively explained by the fact that in order to satisfy the constraint each vertex must be aware of the activity status (x j = 1 or x j = 0) of every other vertex. This implies that the decision variable such as x 1 is related to |F| − 1 other decision variables.
To mitigate the complexity, we split the p-constraint into several constraints X j2S k where fS k g m k¼1 is a partition of the facility set F and all the p k add up to p. Then, the quadratic penalty in (15)  The new objective function is whereQ is a symmetric matrix.

Weights estimation, TOPSIS
The TOPSIS algorithm is widely used in multi-attribute decision-making problems such as supply chain logistics, engineering and facility allocation problems [12,20]. As mentioned before, the algorithm determines the best ranking (weighting) for a set of candidates with certain attributes. The algorithm treats each candidate as a point in the Euclidean space, where the candidate's attributes become spatial coordinates. Each candidate is weighted based on its Euclidean distance to the synthesized ideal solution. The ideal solution is a theoretical candidate with the most desired attributes, and if such a candidate existed, it would get the highest possible ranking. In practice, the ideal candidate may not exist, but it could be synthesized based on the available data and criteria determined by a decision-maker. TOPSIS creates the ideal solution and produces the final ranking (weights) for all candidates based on their Euclidean distance proximity to the ideal solution.
In order to estimate the model's weights w j for j = 1, ‥, |F|, we use TOPSIS from the opensource Python library scikit-criteria [33]. First, we create an alternatives matrix that consists of rows of attributes of each facility. The matrix serves as an input to the TOPSIS algorithm. In this study, we consider three attributes to determine the TOPSIS weights w j . Namely, we use 1. Monthly demand data per facility (demand attribute).
2. Transit facility connectedness to other services (connectedness attribute).
3. Number of neighbouring landmarks that could produce additional demand (landmarks attribute). For example, the landmarks could be hospitals, schools, subways, malls etc.
We note that these data are independent of the parameters used in the optimization model. Table 3 is an example of such a matrix for a problem with three facilities. The W d , W c , W l are user-defined importance parameters of each attribute. In this case, we give the highest preference (W d = .45) to the demand attribute of each facility.
We summarize the TOPSIS workflow in the Fig 3. Note that despite the fact that the facility A has the highest demand, the Fig 3 suggests that facility C has the highest weighting as it has better overall characteristics. We compute weights w j by generalizing this method to all facilities.

Model optimization
Taking advantage of QUBO's equivalence to the transverse Ising model [22] we can view the optimization problem as energy minimization of the Ising objective function. Briefly, the transverse Ising model was originally developed to describe quantum phenomena in

PLOS ONE
Transit facility allocation: Hybrid quantum-classical optimization magnetism. The model consists of discrete variables that represent atomic spins that can be in one of two states (+ 1 or −1). It has been proposed that quantum annealing is an effective physical process that can attain the minimum energy state of a spin system [22,34]. In 2007 a Canadian quantum computing company D-Wave Systems released a quantum device that utilized quantum annealing to solve Ising model-related problems. Since then they have released more powerful hardware suitable for large-scale discrete optimization problems. The D-Wave Quantum Processing Unit (QPU) can be viewed as a heuristic that efficiently minimizes Ising objective functions using a physically-realized version of quantum annealing. In this paper, we utilize D-Wave's quantum hardware along with classical computing heuristics to find the best solution for our problem.

Hybrid solver
It has been shown that decomposing QUBO problems into smaller QUBO sub-problems and combining the resulting solutions can be an effective method for obtaining optimal results [35,36].
To find the best solution, we use an open-source dwave-hybrid Python framework [13]. The hybrid solver attempts to solve the entire problem at once while also solving multiple child sub-problems of the original problem. We integrate the classical heuristic methods tabu search and simulated annealing with quantum annealing into a single system. Both tabu search and simulated annealing work on the entire problem, whereas quantum annealing works on multiple child sub-problems of a smaller size where only high energy impact variables are optimized. The high energy impact variables are the variables whose contribution to the objective function is greater than the rest of the variables. All three heuristics run in parallel branches and combine their best solutions during each iteration k. See Fig 4. The process terminates whenever the time or iteration budget is depleted. We summarize the algorithm structure in Algorithm 1. The use of such a hybrid approach is highly beneficial for several reasons. First, we recall that tabu search and simulated annealing are hill-climbing local search heuristics. During the minimization process, local search heuristics tend to get trapped in suboptimal regions of an optimization landscape. To avoid being trapped in local minima, both methods deploy mechanisms of accepting uphill moves (worsening moves) that could lead to better

PLOS ONE
Transit facility allocation: Hybrid quantum-classical optimization regions. Tabu search can accept uphill moves based on its tabu list and fitness characteristics of a candidate move. By contrast, simulated annealing uses a probabilistic rule called the Metropolis-Hastings criterion [37] for accepting uphill moves. While hill-climbing turns out to be a fairly successful approach, it is not the only way of navigating the optimization landscape. Quantum annealing is also capable of hill-climbing, but it can also perform quantum tunnelling [14]. Intuitively, this can be seen as moving towards the global optimum right through an optimization landscape barrier. Such a tunnelling move is useful when the barrier is so tall that hill-climbing methods would need to accept too many suboptimal moves to traverse the barrier. Therefore, quantum tunnelling is a useful addition that aids the exploration of the solution search space.
Algorithm 1 The hybrid algorithm runs simulated annealing (SA), tabu search (TS) and quantum annealing (QA) on multiple sub-problems. x sa run SA on a problem x t Qx starting with x � . 10: x tb run TS on a problem x t Qx starting with x � . 11: x � choose deterministically or stochastically x sa or x tb . 12: Identify a subset X of size m such that for j 2 X, x j has high energy impact. 13: Formulate a sub-problem y t Qy such that y i ¼ x � i is fixed for i = 2 X. 14: x qa run QA on a problem y t Qy where for j 2 X, y j is a variable.

15:
If needed, repeat Steps 12-14 for different X and track x qa . 16: x � choose the best of x qa , x sa , x tb . 17: k k + 1. 18: end while 19: Return x � . 20: End One might wonder why not entirely rely on quantum computation which can do both hill climbing and tunnelling. The answer is fairly nuanced. In short, the current generation of quantum hardware suffers from noise and errors that occur during quantum computation [38]. Without error correction, the accumulation of errors can significantly degrade the quality of the optimization results and even distort the quantum implementation of a problem [39]. Therefore, relying exclusively on quantum computation is still not entirely feasible. The hybrid approach, on the other hand, takes the best of both worlds where the parallel execution branches create necessary redundancy. Should one system fail due to errors or get trapped, the other branch takes over and provides its solution for the next iteration.

Integrality of solutions under continuous relaxation
In this section, we investigate the properties of a matrixQ in (25) and demonstrate that under mild conditions, the optimization problems (21) and (25) naturally yield integral solutions even if the integrality constraints are disregarded. First,Q is a symmetric matrix. Due to the nature of the problem and our specific formulation choices, the matrixQ has relatively large positive diagonal entries. This can be observed by inspecting (2). Based on a numerical investigation ofQ we learn that the matrix is positive definite whenever p jFj > 0:5 and m j � 2, see Fig 5. In other words,Q is positive definite whenever the ratio of retained facilities to a total number of facilities is more than a half and there are at most 2 competitors per facility. In the case whenQ is not positive definite, we can obtain an equivalent representation of the problem with a perturbed positive definite matrix Q 0 ¼Qðl þ �ÞI where λ is the smallest eigenvalue ofQ and � is chosen such that λ + � > 0 [40].
SinceQ is positive definite, it is informative to consider the continuous relaxation of the problem: We relax the problem by allowing each x j to be a continuous variable restricted to the compact interval [0, 1]. The resulting constraint set is a hyper-cube. Hence, if x is a feasible solution, then x 2 [0, 1] |F| . The positive definiteness ofQ implies that the objective function in (26) is convex. Given that we aim to maximize the objective andQ is positive definite, the optimal solution occurs at one of the corners of the hyper-cube. Since the position of each corner

PLOS ONE
corresponds to a binary string of length |F|, the optimal solution x � is guaranteed to be binary, i.e. x � j 2 f0; 1g for j = 1, . . ., |F|. This means, if we allow 0 � x j � 1, then any suitable continuous optimization solver will yield binary solutions. It is worth noting that even though the objective function is convex we do maximization. Hence, there is no guarantee that the problem will be solved in polynomial time.
Using the continuous optimization solver IPOPT (Interior Point OPTimizer) along with the open-source GEKKO optimization framework [31] we obtain a solution to (26) that is identical to the solution to (25) from the Hybrid solver presented in the previous section. Furthermore, we do another solve using the branch and bound solver APOPT (Advanced Process OPTimizer) and obtain exactly the same solution again. Given the consensus between hybrid, continuous and combinatorial solvers we have a reason to believe that the obtained solution is globally optimal. Moreover, this demonstrates the flexibility of our formulation as it can be solved on a wide variety of discrete and continuous solvers with minimal modifications. In the next section, we investigate the optimality of the solution by presenting an upper bound for the problem.

Upper bound
We would like to ensure that the solutions to the problem are close to a theoretical upper bound for the problem. We claim that for all solutions x 2 {0, 1} |F| that satisfy the partitioned p-constraint (22) the following holds, where λ 1 � . . . � λ p � . . . � λ |F| are eigenvalues of Q defined in (3) and (4) and p is the number of facilities we want to retain. To show this, let |F| = n and ∑ k p k = p. Then We note that the quadratic penalty term is zero because x satisfies the partitioned p-constraint, P j2S k x j ¼ p k for all k. From the definition of Q we know that for all i, Q ii > 0 and Q ij � 0 for i 6 ¼ j. Hence let us drop the off-diagonal terms, Let fQ i k i k Þ g n k¼1 be a non-increasing re-arrangement of fQ ii g n i¼1 . Since the partitioned p-constraint (22) implies p-constraint (2), we have P n i¼1 x i ¼ p. This implies there are n − p decision variables x i that are equal to zero. Hence, Applying the Schur-Horn theorem [41, Theorem 4.3.45] yields Therefore, for all x 2 {0, 1} n satisfying (22), we have This shows that an upper bound can be calculated by summing the p largest diagonal entries or eigenvalues of Q.

Case study
In this section, we apply our framework to the northbound bus route B20 of the transit system in Vancouver, British Columbia. The route connects suburban areas to Vancouver downtown and has 49 bus stops. The following sections present accessibility analysis and route efficiency improvements for transit-rider friendly value p = 40 (20% reduction of bus stops). In addition, we perform a redundancy analysis and show that it is possible to remove up to 40% of the bust stops while maintaining the same accessibility.

Accessibility and efficiency analysis
Using our open-source GIS framework, we create a computer model of the northbound route and surrounding dissemination areas (demand nodes). The model stores information about the population and geographical vector boundaries of each dissemination area, bus stop connectivity to other routes, TOPSIS weights w j , neighbourhoods F j and A j for each bus stop j. As mentioned before, the northbound route B20 has 49 bus stops. We set p = 40, which is roughly 20% reduction of the number of stops. We use a radius of 400 meters for the neighbourhood F j (neighbourhood of competitors around the stop j) and a recommended 500 meters walking distance for the radius of A j (neighbourhood of dissemination areas around the stop j) with α, β = 1. Even though the number of bus stops is reduced by 20%, the retained stops still cover the same number of dissemination areas. This means the recommended walking distance to a bus stop of 500 meters is also retained. Using bus dwell time and boarding/alighting data from Translink we estimate a reduction of travel time by approximately 7 minutes, which is 13% of travel time on the unoptimized route during work hours. Therefore, we reduce the number of bus stops by 20% while preserving the same route accessibility and increasing service effectiveness.

Redundancy analysis
To further analyze the best solutions for different p, we count how many northbound bus stops are within the walking distance of each dissemination area. The counting yields sequences of numbers N p ¼ fn ðpÞ

PLOS ONE
Transit facility allocation: Hybrid quantum-classical optimization manages to find an optimal configuration of retained stops to provide complete coverage. This analysis coincides with the marginal increase of the objective function for p > 30 (Fig 7, yellow line). Further, this suggests that the model effectively handles inter-facility competition.

Additional numerical experiments
This section compares our model solved with the hybrid solver against a SIC type model solved with the APOPT solver. We demonstrate that the proposed methodology yields statistically better optimization results in terms of demand coverage for different numbers of retained facilities.
For this experiment, we randomly generate 30 synthetic problems with 50 facilities and 135 demand nodes. Each facility j has uniformly random weights w j 2 (0, 1]. The synthetic facility and demand node placements are based on the randomly perturbed coordinates of the facilities and demand nodes of the case study route B20. We choose the maximum perturbation offset of 150 meters because it allows generating diverse synthetic problems while still representing a relatively realistic scenario. Each generated instance is solved using two approaches. The first approach A 1 uses the proposed mathematical formulation in (16) and the Hybrid solver presented in 3.5. We set one second time limit on the hybrid solver for each problem instance. The second approach A 2

PLOS ONE
Transit facility allocation: Hybrid quantum-classical optimization uses the SIC type formulation with the APOPT solver. The SIC formulation is as follows ; 1g for all j:

PLOS ONE
In this problem, D denotes the set of all demand nodes, and N i denotes the set of facilities in the neighbourhood of the demand node i. The rest of the notation is defined as in Table 2. The SIC problem in (33) always has a solution for all values of p 2 {1, . . ., |F|} and any location configuration of facilities and demand nodes. We set no time limit on the APOPT solver for each problem instance and allow 100000 maximum iterations to solve the problem. The convergence to a solution happens before the depletion of the iteration budget. We test two hypotheses: the null hypothesis (H 0 ), which states the first approach A 1 covers less demand than A 2 , and an alternative hypothesis (H a ) which states that A 1 covers more demand than A 2 . The demand coverage difference d is the difference between the number of demand nodes covered using A 1 and the number of demand nodes covered using A 2 . If d > 0, then A 1 is better. Using the Wilcoxon signed-rank test [42] and the synthetic dataset, we summarize the results in Table 4.

Comparison with related models
This section qualitatively compares and summarizes the proposed model with the related models reviewed in Section 1.2. Specifically, we compare the proposed model with LSCP [15], p-median [16], MCLP and SIC models [5].
We commence the comparison by emphasizing the fact that the proposed model is primarily designed for quantum optimization. This means the proposed model is the only model amenable for quantum optimization and capable of taking advantage of quantum phenomena such as quantum superposition, quantum tunnelling and entanglement. Another unique differentiating feature of the proposed model is the ability to produce naturally integral solutions even if the integrality constraints are disregarded, and continuous optimization methods are applied. This is not the case for all other related models, as they require heuristic techniques that enforce the integrality of solutions. Finally, it is possible to compute a tight upper bound on the optimal objective function value. This computation can be done efficiently as it is only necessary to compute the sum of p largest diagonal entries of the matrix Q-for further details, see Section 5.
We summarize the comparison of the proposed model with the related models in Table 5.

Results and discussion
Unlike most literature on facility allocation and route optimization, we explore non-linear, physics-inspired mathematical modelling. This allows us to attain an expressive model that captures the most substantial modelling aspects of P-median, Maximal Covering Location and Spatial Interaction Coverage models. During the model derivation process, we demonstrate that it captures complex interactions between facilities and demand nodes, inter-facility competition and distance decay. We show that GIS and demand data can be flexibly integrated using the decision-making analysis techniques such as TOPSIS. To demonstrate our model's computational robustness, we implement a classical-quantum Hybrid solver that harnesses the quantum effects such as superposition and quantum tunnelling while providing redundancy with classical solvers. As a case study, we apply our optimization pipeline on the B20 bus route located in one of the most populated cities of Canada-Vancouver. We demonstrate that the proposed approach reduces the number of bus stop facilities by up to 40% while maintaining the same service coverage. Using the derived theoretical upper bound, we analyze the quality of solutions. The smallest ratio of the best solution to the upper bound is 5954/6291 � 0.95,  [43] suggests that traditional computer devices are struggling to support the growing computational requirements. Due to this, research in alternative computational methods has steadily increased. Quantum computation holds a promise to speed up some of the most demanding optimization processes. This study clearly shows that the proposed optimization framework allows facility planners to consider emerging alternatives to classical computation by experimenting with novel specialized hardware that is capable of harnessing quantum effects.
While the ability to exploit novel quantum phenomena is one of the motivating factors for developing the proposed framework, the most crucial criterion is retaining important modelling aspects present in the predecessor models such as MCLP and SIC. One of the improvements present in our model is the possibility of introducing convexity into the objective function. We argue that whenever the objective function is convex, solutions to the problem are always binary, even if the integrality constraints are disregarded. This allows the proposed model to be treated from the perspective of continuous optimization that is rich in theory and optimization methods.
GIS and demand data can be large and complex [44]. Extracting useful information and deciding what aspects of the data should contribute to the decision-making process might be a great challenge [44,45]. We address this by introducing decision-making analysis which allows facility planners to rank important factors in a formal analytical manner. Using TOPSIS, multiple competing decision criteria can be readily integrated into the optimization model for further analysis.
Although successful in many aspects, our approach suffers from the following limitations: • As the number of retained facilities decreases, there is an inevitable loss in coverage for some demand nodes. When the demand coverage is critically important, and the solutions that lose coverage are deemed infeasible, our model fails as it allows for loss of coverage.
• Due to QUBO formalism, the addition of inequality constraints results in extra auxiliary binary variables. Possibly some inequality constraints may impact an optimization landscape's ruggedness, challenging the hill-climbing heuristics.
• Introducing integer variables that are not binary is a challenge. While not impossible, the result is likely to be inefficient.
We now discuss these limitations. While the loss of coverage may be considered a limitation, the ability to lose coverage guarantees the existence of a solution for all values of p 2 {1, 2, . . ., |F|} and any geographical configuration of facilities and demand nodes. It may be impossible to reduce the number of facilities in many scenarios without losing some demand coverage.
Additional inequality constraints require auxiliary binary variables and need to be incorporated into an objective function. While this may pose a challenge, the non-linear nature of the model allows expressing many complex relations in an alternative fashion; for various methods see [46,47]. Nevertheless, complex inequality constraints may pose mathematical and optimization issues.
The use of integer variables is another challenge from the quantum hardware perspective. The current generation of quantum hardware does not have enough resources to solve non-binary problems. However, recent developments [48] address this issue through various software-based methods.

Conclusion and future works
This paper presents a new framework for optimizing facility allocation problems. The framework aims to reduce the number of facilities while maximizing accessibility. It considers several factors such as geographical Census Program data, facility connectivity, demand data, inter-facility competition and proximity to public landmarks. We augment combinatorial optimization with multi-criterion decision-making analysis to establish data-driven rankings of each facility. This allows for building a parameterizable model that can be easily adapted for different usage scenarios. The proposed mathematical model possesses unique characteristics. It has an efficient computable tight bound that can be used to analyze the quality of solutions. Under mild conditions, problem solutions are naturally integral. The model's equivalence to the Ising spin glass model allows harnessing quantum effects such as superposition and quantum tunnelling.
In future work, we plan to adapt the developed methods to multi-objective problems with distinct objectives capturing various optimization criteria. Since the proposed formalism allows algorithm implementation on fully programmable gate-based quantum computers, we plan to implement a hybrid optimization technique based on Variation Quantum Eigensolver [27] suitable for near-term quantum applications.